Fitness versus Longevity in Age-Structured Population Dynamics 
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We examine the dynamics of an age-structured population model in which the life expectancy 
of an offspring may be mutated with respect to that of the parent. While the total population 
of the system always reaches a steady state, the fitness and age characteristics exhibit counter- 
intuitive behavior as a function of the mutational bias. By analytical and numerical study of the 
underlying rate equations, we show that if deleterious mutations are favored, the average fitness of 
the population reaches a steady state, while the average population age is a decreasing function of 
the overall fitness. When advantageous mutations are favored, the average population fitness grows 
linearly with time t, while the average age is independent of fitness. For no mutational bias, the 
average fitness grows as t 2 ^ 3 . 



I. INTRODUCTION 

The goal of this paper is to understand the role of mu- 
tations on the evolution of fitness and age characteristics 
of individuals in a simple age-structured population dy- 
namics model jl| . While there are many classical models 
to describe single-species population dynamics [^|j3| , con- 
sideration of age-dependent characteristics is a more re- 
cent development J|^-0|. Typically, age characteristics 
of a population are determined by studying rate equa- 
tions which include age-dependent birth and death rates. 
Here we will study an extension of age-structured popula- 
tion dynamics in which the characteristics of an offspring 
are mutated with respect to its parent. In particular, an 
offspring may be more "fit" or less fit than its parent, 
and this may be reflected in attributes such as its birth 
rate and/or its life expectancy. 

In our model, we characterize the fitness of an indi- 
vidual by a single heritable trait - the life expectancy 
n - which is defined as the average life span of an in- 
dividual in the absence of competition. This provides a 
useful fitness measure, as a longer-lived individual has 
a higher chance of producing more offspring throughout 
its life span. We allow for either deleterious or advanta- 
geous mutations, where the offspring fitness is less than 
or greater than that of the parent, respectively (Fig. [I]). 
This leads to three different behaviors which depend on 
the ratio between these two mutation rates. When ad- 
vantageous mutation is favored, the fitness distribution of 
the population approaches a Gaussian, with the average 
fitness growing linearly with time t and the width of the 
distribution increasing as t 1 / 2 . Conversely, when delete- 
rious mutation is more likely, a steady-state fitness distri- 
bution is approached, with the rate of approach varying 
as t~ 2 / 3 . When there is no mutational bias, the fitness 
distribution again approaches a Gaussian, but with av- 
erage fitness growing as t 2 ! 3 and the width of the distri- 
bution again growing as t 1 / 2 . 

In all three cases, the average population age reaches 
a steady state which, surprisingly, is a decreasing func- 
tion of the average fitness. Thus within our model, a 
more fit population does not lead to an increased individ- 
ual lifetime. Qualitatively, as individuals become more 



fit, competition plays a more prominent role and is the 
primary mechanism that leads to premature mortality. 
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FIG. I. Schematic illustration of the model. An individual 
with fitness, or intrinsic life expectancy n, continuously ages 
(horizontal arrow). The heavy dots signify individual birth 
events. At each birth event, an offspring (of age zero) is pro- 
duced, whose intrinsic lifetime is either n+1, n, or n-1, with 
relative rates b+, b, and 6_, respectively. The individual dies 
either by aging or by competition ( x ) . 

In the following two sections, we formally define the 
model and outline qualitative features of the population 
dynamics. In Sees. IV- VI, we analyze the three cases of 
deleterious, advantageous, and neutral mutational biases 
in detail. We conclude in Sec. VII. Various calculational 
details are provided in the Appendices. 



II. THE MODEL 

Our model is a simple extension of logistic dynamics 
in which a population with overall density N(t) evolves 
both by birth at rate b and death at rate 7 AT. Such a 
system is described by the rate equation 



N = bN — -1N 2 



(1) 



with steady-state solution = 6/7. Our age- 

structured mutation model incorporates the following ad- 
ditional features: 

1. Each individual is endowed with given life ex- 
pectancy n. This means that an individual has a 
rate of dying which equals \/n. 
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2. Death by aging occurs at a rate inversely propor- 
tional to the life expectancy. 

3. Individuals give birth at a constant rate during 
their lifetimes. 

4. In each birth event, the life expectancy of the new- 
born may be equal to that of its parent, or the life 
expectancy may be increased by or decreased by 1. 
The relative rates of these events are b, b+, and 
respectively. 

Each of these features represent idealizations. Most 
prominently, it would be desirable to incorporate a re- 
alistic mortality rate which is an increasing function of 
age However, we argue in Sec. VII that our basic 

conclusions continue to be valid for systems with realistic 
mortality rates. 

To describe this dynamics mathematically, we study 
C n (a,t), the density of individuals with life expectancy 
n > 1 and age a at time t. We also introduce P n (t) — 
J C n (a,t) da, which is the density of individuals with 
life expectancy n and any age at time t. Finally, the to- 
tal population density is the integral of the population 
density over all ages and life expectancies, 



N{t) 



oo poo °° 

/ C n (a,t)da=J2 p n(t)' 

7i=l ^° n=l 



(2) 



According to our model, the rate equation for C n (a,t) 



jN(t) + -)C n (a,t). (3) 



The derivative with respect to a on the left-hand side ac- 
counts for the continuous aging of the population ||§). 
On the right-hand side, ~fNC n is the death rate due to 
competition, which is assumed to be independent of an 
individual's age and fitness. As discussed above, the mor- 
tality rate is taken as age independent, and the form 
C n /n guarantees that the life expectancy in the absence 
of competition equals n. Because birth creates individ- 
uals of age a — 0, the population of newborns provides 
the following boundary condition for C n (0,t), 



= 16- 7-ZV — — ) P n + b+Pn-l + b-Pn+1- (5) 



dP n 
dt 



This describes a random-walk-like process with state- 
dependent hopping rates in the one-dimensional fitness 
space n. Notice the hidden non-linearity embodied by the 
term -fNP n , since the total population density N(t) — 
J2n Pn{t)- From Eq. (||), we find that N(t) obeys a gen- 
eralized logistic equation 



b--jN)N- V — -6_Pi. (6) 
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C n (0,t) - bP n (t) + 6+P„_i(t) + 6-P„+i(i). 



(4) 



FIG. 2. Illustration of the random walk in fitness space that 
underlies the behavior of P n . The term (b — ~/N — —)P n < 
represents a state-dependent population loss (dashed arrows) 
which decreases for larger n. 

The three different dynamical regimes outlined in the 
introduction are characterized by the relative magnitudes 
of the mutation rates 6+ and b- . The main features of 
these regimes are: 

• Subcritical case. Here b + < 6_, or deleterious 
mutations prevail. The fitness of the population 
eventually reaches a steady state. 

• Critical case. Here b + = &_, or no mutational 
bias. The average fitness of the population grows as 
i 2 / 3 and the width of the fitness distribution grows 
as t 1 ' 2 . 

• Supercritical case. Here 6+ > 6_, or advanta- 
geous mutations are favorable. The average fitness 
grows linearly in time and the width of the fitness 
distribution still grows as t 1 / 2 . 

In all three cases, the total population density N and 
the average age A, defined by 



Finally, the condition Pq = follows from the require- 
ment that offspring with zero life expectancy cannot be 
born. 



III. BASIC POPULATION CHARACTERISTICS 

Let us first study the fitness characteristics of the pop- 
ulation and disregard the age structure. The rate equa- 
tion for P n {t) is found by integrating Eq. (|^) over all ages 
and then using the boundary condition Eq. (||) to give 



1 

A= N^2 I daaC n (a) 



(7) 



saturate to finite values. The steady state values of N 
and A may be determined by balance between the total 
birth rate B = b + b + + 6_ and the death rate 7^V due to 
overcrowding. For example, in the critical and supercrit- 
ical cases, there are essentially no individuals with small 
fitness, so that the last two terms in Eq. (g) may be ne- 
glected. Then the steady-state solution to this equation 
is simply 
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N = —. 

7 



(8) 



This statement also expresses the obvious fact that in the 
steady state the total birth rate B must balance the total 
death rate 7 A. (For fit populations, the death rate due 
to aging is negligible.) Similarly, the average age may 
be inferred from the condition it must equal the average 
time between death events. Thus 



A = 



1 

77V 



1 

B' 



(9) 



The behavior of the average age in the subcritical case is 
more subtle and we treat this case in detail in the section 
following. 



IV. THE SUBCRITICAL CASE 



When deleterious mutations are favored (6_ > b + ), the 
random-walk-like master equation for P n contains both 
the mutational bias towards the absorbing boundary at 
the origin, as well as an effective positive bias due to the 
1 jn term on the right-hand side of Eq. (||) . The balance 
between these two opposite biases leads to a stationary 
state whose solution is found by setting P n = in Eq. (||). 
To obtain this steady state solution, it is convenient to 
introduce the generating function 



straightforwardly by a direct integration of Eq. (|T^) , only 
one member of this family is the correct one. To deter- 
mine this true solution, we must invoke additional argu- 
ments about the physically realizable value of N for a 
given initial condition. 

An upper bound for N may be found from the steady- 
state version of Eq. (||), 



(B - 7 A)A = V — + b-P ± 



(14) 



Since the right-hand side must be non-negative, this pro- 
vides the bound 7A < B. On the other hand, we may 
obtain a lower bound for N by considering the master 
equation for P n in the steady state. For n — ► oo, we 
may neglect the P n /n term in Eq. (^) and then solve the 
resulting equation to give P n = A + A™ + A_A™ , where 



A± = \jN -b± y/(l N - b) 2 - 46+b_l /2b 



(15) 



For P n to remain positive, X± should be real. Hence we 
require 7 A > b + 2 yjb + b. 
N must lie in the range 



We therefore conclude that 



b + 2 v / b+b- < 7 A < B. 



(16) 



F(x)=J2 P nX n - 1 . 



(10) 



Multiplying Eq. (Q) by x" 1 and summing over n gives 



= (b-jN)F- 



b+xF + b- 



n=l 



-Pi 



The term involving P n /n is simplified by using 

n _n— 1 



— a 

n 



1 r 

- / F(y)dy. 
x Jo 



(11) 



(12) 



Multiplying Eq. (|Tl|) by x and differentiating with respect 
to x gives 



F'(x) _ jN-b+1- 2b + x 
F(x) ~ b+x 2 - (jN - b)x + b- ' 



(13) 



where the prime denotes differentiation. 

As in the case of the master equation for P n , this dif- 
ferential equation for F has a hidden indeterminacy, as 
the total population density N — F(x — 1) appears on 
the right-hand side. Thus an integration of Eq. (13), 
subject to the boundary condition F(l) = N, actually 
gives a family of solutions which are parameterized by the 
value of N. While the family of solutions can be obtained 



While the foregoing suggests that N lies within a finite 
range, we find numerically that the minimal solution, 
which satisfies the lower bound of Eq. ([l6|), is the one 
that is generally realized (Fig. ||). This selection phe- 
nomenon is reminiscent of the corresponding behavior 
in the Fisher-Kolmogorov equation and related reaction- 
diffusion systems 0] , where only the extremal solution is 
selected from a continuous range of a priori solutions. 



1.2 
1.0 
0.8 
£ 0.6 
0.4 
0.2 
0.0 



0.0 



0.5 



1.0 



1.5 



2.0 



FIG. 3. Minimal steady sta te valu e of the total density N 
versus mutational bias n = yJb+JbZ. Here 7=1,6 = 0, and 
b+ + b- = 1, so that the total birth rate B = b + b+ + 6_ is 
fixed. For n > 1, N sticks at the value of unity. 
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To understand the nature of this extremal solution in 
the present context, notice that with the bounds on N 
given in Eq. ([l6|), A + lies within the range [//, 1), where 



m= i r 



(17) 



is the fundamental parameter which characterizes the 
mutational bias. Consequently the steady-state fitness 
distribution decays exponentially with n, namely P n ~ 
A" . When the total population density attains the min- 
imal value -ZV min = (b + 2y / b + b-)/j, A + also achieves 
its minimum possible value A" lln = /x, so that the fit- 
ness distribution has the most rapid decay in n for 
the minimal solution. Based on the analogy with the 
Fisher-Kolmogorov equation |l , we infer that there are 
two distinct steady-state behaviors for P n as a function 
of the initial condition P n (0). For any P n (0) with ei- 
ther a finite support in n or decaying at least as fast 
as fi n , the extremal solution P n ~ /i™ is approached 
as t — > co. Conversely, for initial conditions in which 
P„(0) decays more slowly than zx n , for example as a", 
with a in the range (zx,l), the asymptotic solution also 
decays as a™. Correspondingly, Eq. (|J) in the steady 
state predicts a larger than minimal population density 
N = (b + b-a + b+a.- 1 )/^. 
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FIG. 4. N(t) versus t 2 ^ z in the subcritical case for 6 = 0, 
b + — 1, 6_ = 2, and 7 = 0.5 for the initial condi- 
tions: (i) P„(t = 0) = 0.1 for 1 < n < 10 (o), (ii) 
P n (0) = a", with a = (1 + fi)/2 (A), and (iii) P n (0) = 1/n 2 
for (V). Asymptotically, the data for N(t) approach the 
respective theoretical values of iV m i n = 4\/2 ~ 5.6568, 
AT(co) = (2a + a _1 )/7 w 5.7574, and iV(oo) = B/7 = 6. 
The rate of approach is t~ 2 ^ 3 in the first case and faster than 
t~ 2 ^ 3 in the latter two cases. These calculations use the full 
machine precision of 10 -308 . 

We also find that the extremal and the non-extremal 
solutions exhibit different relaxations to the steady state. 
For those initial conditions which evolve to the extremal 
solution, the deviation of N and indeed each of the P n 



from their steady state values decay as t~ 2 ^ 3 , while for all 
other initial conditions, the relaxation to the steady state 
appears to follow a t~ x power law decay. The power-law 
approach to the steady state is surprising, since the over- 
all density obeys a logistic- like dynamics, N = bN—jN 2 , 
for which the approach to the steady state is exponential. 
These results are illustrated in Fig. [| which shows the 
asymptotic time dependence of N(t) based on a numer- 
ical integration of Eq. (||) with the fourth-order Runge- 
Kutta algorithm ||. The demonstration of the i -2 / 3 re- 
laxation to the extremal solution relies on a correspon- 
dence to the transient behavior in the critical case. This 
is presented in Appendix B. 
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FIG. 5. Behavior of N(t) versus t 2 ^ z for precision equal to 
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10" 



10" 



and 10" 



10" 1UU , 10" 

for the case 6+ = 16, b- = 1, b = 0, and 7 = 0.5. In the limit 
of infinite precision N(t) — + 16 as t — > 00. 

A disconcerting feature of the numerical calculation for 
N(t) is the small disagreement between the numerically 
observed values of the steady-state population density 
and the expected theoretical prediction (Fig. ||) . This dis- 
crepancy arises from the finite computer precision which 
causes very small values of P n to be set to zero. To con- 
firm this, we changed the computer precision from 10 _ 
to the full machine precision of 1(T 308 (Fig. |). As the 
precision is increased, saturates to progressively higher 
values and approaches the theoretical prediction. A sim- 
ilar precisions-dependent phenomenon has been observed 
in the context of traveling Fisher-Kolmogorov wave prop- 
agation 

For the relevant situation where the density Af takes 
the minimal value, we may rewrite Eq. (13) as 



EL 

~F 



2fi 



1 



1 — (ix £>_(1 — (jlx) 2 

Integrating from x — 1 to x and using F(l) = N gives 
2 



(18) 



1 



1 



1 



\ 1 — [ix 1 — /i 



(19) 
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One can now formally determine P n by expanding F(x) 
in a Taylor series. For example, 



Pi = N(l - n) 2 exp 



1 



Mi-/*) 



P 2 =N(l-ny \ 2n+—) <>x]> 



For many applications, however, there is no need to deal 
with these unwieldy expressions. As we now discuss, the 
overall fitness or age characteristics of the population can 
be obtained directly from the generating function with- 
out using the explicit formulae for the P n - 



A. Fitness characteristics 

Consider the average fitness of the population 

< n > = ]vE" p «> (20) 

n=l 

which can be expressed in terms of the generating func- 
tion as 



1 dF 



1. 



x=l 



N dx 

From Eq. ( |l9| ) we thereby obtain the average fitness 
. . 2fi 1 



1-H 6_(1- M ) 2 



1. 



(21) 



(22) 



As one might anticipate, the average fitness diverges as 
fi — ► 1 from below, corresponding to the population be- 
coming mutationally neutral. To determine the disper- 
sion of the fitness distribution we make use of the relation 



(«(«-i)) = ^E n ( n - 1 ) p « = ^ 



1 d 2 {xF) 



N dx 2 



(23) 



Substituting Eqs. ( |l9|) and also Eq. (^2|) then gives 



1 + 



G// 



+ 



3(1 + M ) 



Thus the dispersion a 2 = (n 2 ) — (n) 2 in the fitness dis- 
tribution is 



2// 



+ 



(1-M) 2 6-(l- M ) 3 



(24) 



As the mutational bias vanishes, fj, — ► 1, the average fit- 
ness and the dispersion diverge as (n) ~ &I (1 — A 4 ) 2 
and cr ~ y/2/b-(l — /i)~ 3 / 2 . Thus these two moments 
are related by cr ~ (n) 3 / 4 . As we shall see in Sec. VI, 
this basic relation continues to hold in the critical case. 



B. Age characteristics 

In the steady state, we solve Eq. (|J) to give the con- 
centration of individuals with age a and fitness n 



C n {a) = P n \jN+-J exp 
The average age of the population is 

^ 00 POO 

A = — ^2 / daaC n (a) 

-. oo 

-Y 



jN + — \ a 

n 



(25) 



P, 



N ^ jN + n- 1 ' 



(26) 



where the second line is obtained by using Eq. (|25|) . This 
expression can be rewritten directly in terms of the gen- 
erating function by first noticing that 

pi pi 00 

/ x v F{x)dx = / V P„x" +l/ - 1 dx 
Jo Jo „_, 



E 



Pn 



n + v 



(27) 



Thus we re-express Eq. (g6|) in a form which allows us 
to exploit Eq. (|27|). After several elementary steps, we 
obtain 



.4 



111 



E 



Pn 



jN N (7A) 2 ^ n + (7A)- 1 
111/" 1 

dx x~ F{x). 



7 A N (7 AO 2 J 



(28) 
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FIG. 6. Average age A of t he stea dy state population ver- 
sus mutational bias (i — yb+ /b- . The total birth rate 
B = b + b+ + b- is held fixed throughout, with b — 0, 
6+ + b- = 1, and 7 given by the extremal steady-state so- 
lution 77V = b + 2^/b+b- . For fj, > 1, the average age sticks 
at the value of unity, while for fj, — > 0, A 2. 
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This expression should be compared with the result 
for the critical and supercritical cases, namely A — 
(jN)^ 1 = B^ 1 (see Eq. (||)). In the subcritical case, 
77V < B and the above two expressions A m [ n — B^ 1 and 
Anax = (7^0 1 provide lower and upper bounds for the 
average age. This is proved in Appendix A. Fig. ^| shows 
the surprising feature of Eq. (|2^) that the average age de- 
creases as the population gets fitter! We also see that the 
average age of the least fit population (fi — > 0) is twice 
that of the increasingly fit populations in the critical and 
supercritical cases. We now demonstrate this fact. To 
provide a fair comparison we take the total birth rate 
rate to be equal to unity in both cases and also choose 
b = for simplicity. For fit populations (critical and su- 
percritical cases), the average age is simply A = B^ 1 = 1. 
For the least fit population fi — > 0, and correspondingly 
N — > 0. In this limit, we may write Eq. ( p6| ) as, 

00 00 

71=1 ' 7i=l 

On the other hand, from Eq. (|22| ) we have 

In) « 1 + -r- « 2. (30) 

o_ 

The relation A = (n) is natural for the least fit popula- 
tion, as the total density is small and competition among 
individuals plays an insignificant role. Thus the average 
age may be found by merely averaging the intrinsic life 
expectancy of the population. Intriguingly, in this limit 
the average individual in the least fit population lives 
twice as long as individuals in relatively fit populations. 

It is also worth noting that in the limit of a minimally 
fit population (u — > 0) we can expand the generating 
function in Eq. (19) systematically. We thereby find 



that the density P n exhibits a super-exponential decay, 
P n = JVe-V(n - 1)!. 



V. THE SUPERCRITICAL CASE 

When advantageous mutations are favored, the mas- 
ter equation for P n , Eq. (||), can be viewed as a random 
walk with a bias towards increasing n. Because there is 
no mechanism to counterbalance this bias, the average 
fitness grows without bound and no steady state exists. 
As in the case of a uniformly biased random walk on a 
semi-infinite domain, the distribution of fitness becomes 
relatively localized in fitness space, with the peak drift- 
ing towards increasing n with a velocity V — b + — &_. 
Since the fitness profile is non-zero only for large n in 
the long time limit, it is appropriate to adopt continuum 
approach. We therefore treat n as continuous, and derive 
the continuum limits of Eqs. (|J) and (||). For the time 
evolution of the fitness distribution P(n, t), we obtain the 
equation of motion 



d_ 

dt 



v 



ch 



P = [B-jN- 



P + D 



_ d 2 P 
dn 2 



(31) 



This is just a convection-diffusion equation, supple- 
mented by a birth/death term. Here the difference be- 
tween advantageous and deleterious mutations provides 
the drift velocity V = b + — &_ , and the average mutation 
rate D = (b + + &_)/2 plays the role of diffusion constant. 
For the total population density we obtain 



dN 
~dt 



= (B - 7 AT) N - [ 
Jo 



dn 



P(n,t) 



(32) 



To determine the asymptotic behavior of these equa- 
tions, we use the fact that the fitness distribution be- 
comes strongly localized about a value of n which in- 
creases as Vt. Thus we replace the integral in Eq. ( |3^ ) 
by its value at the peak of the distribution, N/Vt. With 
this crude approximation, Eq. (p3) becomes 



™ = ( B -~ / N-±)N. 
dt V Vt 



(33) 



Thus we conclude that the density approaches its steady 
state value as 



7A — > B 



1 

Vt' 



(34) 

as well as the rate of 



This provides both a proof Eq. 
convergence to the steady state. 

We now substitute this asymptotic behavior for the 
total population density into Eq. (J3l]) to obtain 



(35) 



at on I \Vt n J on z 



Notice that the birth/death term on the right hand side 
is negative (positive) for subpopulations which are less 
(more) fit than average fitness Vt. This birth/death term 
must also be zero, on average, since the total population 
density saturates to a constant value. Moreover, this 
term must be small near the peak of the fitness distri- 
butions where n ~ Vt. Thus as a simple approximation, 
we merely neglect this birth/death term and check the 
validity of this assumption a posteriori. This transforms 
Eq. (pq) into the classical convection-diffusion equation 
whose solution is 



P(n,t) 



N 



exp 



(n - Vtf 



ADt 



(36) 



This basic results implies that the fitness distribution in- 
deed is a localized peak, with average fitness growing 
linearly in time, (n) = Vt, and width growing diffu- 
sively, o = \ / 2Dt. We now check the validity of drop- 
ping the birth/death term in Eq. (|35|). Near the peak, 
|n — V£j ~ VDt, so that the birth/death term is of or- 
der i~ 3 / 2 x P. On the other hand, the other terms in 



G 



Eq. ( |35| ) are of order t -1 x P, thus justifying the neglect 
of birth/death term. 

We now turn to the age characteristics. Asymptoti- 
cally, the density of individuals with given age and fit- 
ness changes slowly with time because the overall density 
reaches a steady state. Consequently, the time variable 
t is slow while the age variable a is fast. Physically this 
contrast reflects the fact that during the lifetime of an 
individual the change in the overall age characteristics of 
the population is small. Thus in the first approximation, 
we retain only the age derivative in Eq. (Q). We also 
ignore the term C n /n, which is small near the peak of 
the asymptotic fitness distribution. Solving the result- 
ing master equation and using the boundary condition of 
Eq. (Eh we obtain 



1 N 2 



VAitDt 



exp 



jNa - 



{n - Vtf 
ADt 



(37) 



Integrating over the fitness variable, we find that the age 
distribution C(a,t) — J dnC n {a,t) has the stationary 
Poisson form 



C(a) 



lN 2 e -lNa_ 



(38) 



From this, the average age is A = (jN)^ 1 — B^ 1 in 
agreement with Eq. (||). As discussed in Sec. Ill B, the 
surprising conclusion is that the average age in the super- 
critical case is always smaller than that in the subcritical 
case. 



VI. THE CRITICAL CASE 

We now consider the critical case where the rates of ad- 
vantageous and deleterious mutations are equal. Without 
the 1/n term and with 7iV — b = b + + Eq. (|J) be- 
comes the master equation for an unbiased random walk 
on the semi-infinite range n > 0. Due to the 1/n term, 
the system has a bias towards increasing n which van- 
ishes as n — > oo (see Fig. |^). Thus we anticipate that the 
average fitness will grow faster than t 1 / 2 and slower than 
t. Hence we can again employ the continuum approach 
to account for the evolution of the P n . In this limit, the 
corresponding master equation for P(n, t) becomes 



dP 



D 



•yN — — 
n 



P + D 



. d 2 P 
dn 2 



(39) 



Numerically, we find (n) ~ t 2 / 3 , while the dispersion still 
grows as t 1 / 2 , that is, as a ~ \ft. Thus these two quanti- 
ties are related by a ~ (n) 3 / 4 , as derived analytically for 
the subcritical case. 

To provide a more quantitative derivation of the above 
scaling laws for (n) and a, as well as to determine the 
fitness distribution itself, we examine the equation for 
P(n,t). First note that the total population density still 



obeys Eq. (|32|), as in the supercritical case. Under the 
assumption that the fitness distribution is relatively nar- 
row compared to its mean position, a result which we 
have verified numerically, we again estimate the integral 
on the right-hand side of Eq. (32) to be of the order of 
N/(n). This leads to 



7 A — > B 



1 

M 



Substituting this into Eq. (p9|) yields 



dP 



P + D 



d 2 P 
dn 2 



(40) 



(41) 



Given that the peak of the distribution is located near 
n w (n), it proves useful to change variables from (n, t) to 
the comoving co-ordinates (y = n — (n) , t) to determine 
how the peak of the distribution spreads. We therefore 
write the derivatives in the comoving coordinates 



d_ 

dt 



d_ 

dt 



d(n) d_ 
dt dy' 



d_ 

dn 



d_ 



and expand the birth/death term in powers of the devi- 
ation y = n — (n) 



1 

(n) 



1 

n 



y 



y 



Now Eq. §n\) becomes 

dP d(n) dP _ 
dt dt dy ( 



P 



P + D 



d 2 P 



dy 2 



(42) 



Let us first assume that the average fitness grows faster 
than diffusively, that is, (n) ^> \[t. With this assump- 
tion, the dominant terms in Eq. (M2) are 



d(n) dP 
dt dy 



(43) 



Using 



These terms balance when (n)(ty)~ 
this scaling and balancing the remaining sub-dominant 
terms in Eq. gives y ~ \ft. The combination of 
these results yields (n) ~ t 2 / 3 . This justifies our initial 
assumption that (n) S> Now we write (n) = (ut) 2 ^ 3 , 
with u of order unity, to simplify Eq. (W3J) to 



dP 
dy 



32/ 



P. 



In terms of n 



P(n,t) = N 



2u 2 t 

n) the solution is 



4nu 2 t 



exp 



3 [n - {ut) 2 ' 3 } 
A^t 



(44) 



(45) 



Thus the fitness distribution is again Gaussian, as in the 
supercritical case, but with the average fitness growing 
as (n) = {ut) 2 / 3 . Finally, we determine u = \/3D by 
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substituting (n) = (ut) 2 ^ 3 in Eq. ( [42| ) and balancing the 
sub-dominant terms. 

The age distribution in the critical case can be ob- 
tained in similar manner as in the supercritical case. The 
approximations that were invoked to determine the age 
distribution in the supercritical case still apply. Conse- 
quently, the asymptotic form for C„(a) is still given by 
Eq. (|37|), and this gives the same expression for C(a) af- 
ter integrating over n, as in Eq. (|38|) . Hence the average 
age is again B -1 , as in Eq. (^|). 

VII. SUMMARY AND DISCUSSION 

We have introduced an age-structured logistic- like pop- 
ulation dynamics model, which is augmented by fitness 
mutation of offspring with respect to their parents. Here 
fitness is quantified by the life expectancy n of an individ- 
ual. We found unusual age and fitness evolution in which 
the overall mutational bias leads to three distinct regimes 
of behavior. Specifically, when deleterious mutations are 
more likely, the fitness distribution of the population ap- 
proaches a steady state which is an exponentially decay- 
ing function of fitness. When advantageous mutations are 
favored or when there is no mutational bias, a Gaussian 
fitness distribution arises, in which the average fitness 
grows as (n) = Vt and as (n) = respectively. 

Paradoxically, the average age of the population is 
maximal for a completely unfit population. Conversely, 
individuals are less long-lived for either positive or no mu- 
tational bias, even though the average fitness increases 
indefinitely with time. That is, a continuous "rat-race" 
towards increased fitness leads to a decrease in the aver- 
age life span. As individuals become fit, increased com- 
petition results in their demise well before their intrinsic 
lifetimes are reached. Thus within our model, an increase 
in the average fitness is not a feature which promotes 
longevity. 

Our basic conclusions should continue to hold for the 
more realistic situation where the mortality rate increases 
with age [HQJ§]. The crucial point is that old age is 
unattainable within our model, even if individuals are 
infinitely fit. When the mutational bias is non-negative, 
old age is unattainable due to keen competition among 



fit individuals, while if deleterious mutations are favored, 
age is limited by death due to natural mortality. In ei- 
ther case, there are stringent limits on the life expectancy 
of any individual. To include an age-dependent mortal- 
ity into our model, we may write the mortality term 
f n (a)C n (a,t) instead of n~ 1 C n {a,t) in Eq. (||), where 
/ n (a) is the mortality rate for individuals of age a. Sim- 
ilarly, the term n~ 1 P n in Eq. (|^) should be replaced by 
J daf n (a)C n (a,t). However, these generalized terms play 
no role for large n, since f n (a) is a decreasing function of 
n and old age is unattainable. 
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APPENDIX A: BOUNDS FOR THE AVERAGE AGE 

The upper bound, A < (jN) -1 , follows immediately from Eq. (p8|), so we just prove A > A m i n . We have 
A -B- 1 1 



6 + Ml + /i 2 ) 

A=- — — - / cfa^H^bTTT I _ ZJL ) exp 



6 + 26_/i (b + 2fo_/i) 2 J \1 — fix J \ b-fi \1 — fix l — fi / 

Using these expressions and performing elementary transformations we reduce the inequality A > A min to 
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dx 



Let us now introduce the variable 



exp 



1/1 L_H & + 26_/i 

b-fJ,\l-fix 1-fijj 6 + 6_(l+/x 2 ) 



1 



1 



so that dv = dx/b-(l — /.ix) 2 , which varies in the range [0, V], with V 

1- V- Y v 



dve~ v 



1 

b-[i \ 1 — [ix 1 — /! / 

g ■ This simplifies Eq. (Al) to 

6 + 26_^ 



(Al) 



(A2) 



We now use the inequality 



1 - V-^nv 



1-p 



< 



6 + 6_(l + /z 2 )' 



1-9 



< e 



(A3) 



(A4) 



which holds for < q < p < 1 and v > 0. This inequality is readily proven by taking the logarithm on both sides and 
using the expansion ln(l — u) = — ^2k>i U V^- Now we apply Eq. (A4) to the integrand in (A3) and then replace the 
upper limit V in the integral by oo to give 



dv e 



1 - V^v 



< J dv exp <^— v — v 

,oo r 

< / exp < — v 



6 + 26_/i 

6 + 6_(i + ^ 2 ) 



6 + 26_^ 



b + 2b_fi 
6 + 6_(l + M 2 )' 



This completes the proof. 

The lower bound A m i n turns out to be very accurate in the case when mutations are slightly deleterious. To see 
this let us write 6+ = 1, b- = (1 + e) 2 , where e<l. Replacing x by the transformed variable v — (T 1 — (1 + e — a;)" 1 
recasts the integral Eq. (p2q) as 



dv e~ v 1 



1 — ev 



(A5) 



We now expand the integrand, 



1 — ev 



6 + 2 + 2e 6 + 2 + 2e 



0(e 4 



replace the upper limit in the integral Eq. ( |A5| ) by oo, and compute the resulting simple integrals explicitly to obtain 
a series expansion in e for the average age. Together with analogous expansions for A max and A we have 

1 



6 + 2 + 2e 



A ■ — A 
A = A n 



1 hO(e 6 ) 

(6 + 2 + 2e) 2 + (6 + 2 + 2e) 3 + ( 1 



2e 5 



(6 + 2 + 2e) 2 (6 + 2 + 2e) 3 (6 + 2 + 2e) 3 
Thus the difference between the exact value and A m - ln is of order e 5 . 



0(e 6 



(A6) 
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APPENDIX B: TRANSIENT BEHAVIOR OF THE TOTAL DENSITY 



Numerically, we found that in the subcritical case the total population density approaches the steady state value 
Noo — (b + 2y / b + b-)/j from below with a deviation that vanishes as t~ 2 ^ 3 . We now explain this behavior by 
constructing a mapping between this transient behavior in the subcritical case and the transient behavior in the critical 
case. We start with the basic rate equation, Eq. (||). We may remove the term "fNP n through the transformation 



which simplifies Eq. (|J) to 



dt 



l (t)=P n (t)exp\ 1 / dt'N(t') 



= b-- 



b + Q n -i + b-Q n+1 . 



(Bl) 



(B2) 



Next, the steady state behavior P n ~ /i™ suggests replacing the Q n by R n (t) = ^ n Q n {t). This also removes the 
asymmetry in the birth terms and gives 



dRn 

dt 



b Rn + 6*(i?ri-l + Rn+l) 

n ' 



(B3) 



where we use the shorthand notation b* = ^Jb+bZ. 

One cannot use the continuum approximation to determine the steady-state solutions for P n or Q n . However, the 
continuum approximation is appropriate for the R n . Then Eq. (B3) reduces to 



dR 
dt 



= ( 6 + 26* — 



R + h 



8 2 R 
dn 2 ' 



(B4) 



which is very similar to Eq. (pill). Hence we expect that the distribution of R n is peaked around (n) ~ (3fo* ) 1 / 3 ^ 2 / 3 . 
It proves convenient to make this scaling manifest. To this end we change variables once more, 



*„(/> = /?„(/) <'Xj><| -(b 2/., !/ + ( y 



(B5) 



to get 



dS_ 
dt 



1 1 

(n) n 



S + b : 



d 2 s 

dn 2 



Repeating the procedure of Sec. V we determine the asymptotic solution to Eq. (pq) as 



S n (t) 



1 



exp 



4M 



(B6) 



(B7) 



To find the asymptotics of the total population density let us compute ^ Q n {t)- First, (Bl) can be expressed as 



°° (ft 
V Q n (t) = N(t) exp 7 dt' N(t') 

n=l L Jo 



On the other hand, 



Qn(t) = V nR n{t) = exp <^ (b + 2K)t 



n=l 



to 



1/3- 



(B8) 



(B9) 



In the last s um , the factor fj, n suggests that only terms with small n contribute significantly. Although the asymptotic 
expression (B7) is formally justified only in the scaling region, where \n — (n)\ ~ y/b*t, the continuum approach 
typically provides a qualitatively correct description even outside this region. Therefore we take Eq. (B7) to estimate 
S n for small n. We find 
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£M n Sn(*)~exp|-c(Q 1/8 | 



(BIO) 



where we use (n) ~ i 2 / 3 , as in the critical case, and C is a constant. By substituting Eq. ( BIO ) into Eq. ( |B9|) we 
obtain 



Q„(t) - exp J (6 + 26,)t - (1 + C) 



n=l 



9f 



1/3' 



(Bll) 



Combining Eqs. (B8) and (Bll) we arrive at the asymptotic expansion 



7 / dt'N(t') = (b + 2b*)t-{l + C) - + 



9f 



1/3 



(B12) 



which implies 



7-/V(i) = b + 2b* - const x t~ 2 / 3 + . . 



(B13) 
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